Symmetric and asymmetric solitons in linearly coupled Bose-Einstein 
condensates trapped in optical lattices 



Arthur Gubcskys and Boris A. Malomcd 
Department of Interdisciplinary Studies, 
School of Electrical Engineering 
Faculty of Engineering, Tel Aviv University 
Tel Aviv 69978, Israel 

We study spontaneous symmetry breaking in a system of two parallel quasi-one-dimensional traps 
(cores), equipped with optical lattices (OLs) and filled with a Bose-Einstein condensate (BEC). The 
cores are linearly coupled by tunneling (the model may also be interpreted in terms of spatial 
solitons in parallel planar optical waveguides with a periodic modulation of the refractive index). 
Analysis of the corresponding system of linearly coupled Gross-Pitaevskii equations (GPEs) reveals 
that spectral bandgaps of the single GPE split into subgaps. Symmetry breaking in two-component 
BEC solitons is studied in cases of the attractive (AA) and repulsive (RR) nonlinearity in both traps; 
the mixed situation, with repulsion in one trap and attraction in the other (RA), is considered too. 
In all the cases, stable asymmetric solitons are found, bifurcating from symmetric or antisymmetric 
ones (and destabilizing them), in the AA and RR systems, respectively. In either case, bi-stability is 
predicted, with a nonbifurcating stable branch, either antisymmetric or symmetric, coexisting with 
asymmetric ones. Solitons destabilized by the bifurcation tend to rearrange themselves into their 
stable asymmetric counterparts. In addition to the fundamental solitons, branches of twisted (odd) 
solitons in the AA system, and twisted bound states of fundamental solitons in both AA and RR 
systems, are found too. The impact of a phase mismatch, A, between the OLs in the two cores is 
also studied. It is concluded that A = n/2 only mildly deforms the picture, while A = it changes 
it drastically, replacing the symmetry-breaking bifurcations by pseudo-bifurcations, with the branch 
of asymmetric solutions asymptotically approaching its symmetric or antisymmetric counterpart 
(in the AA and RR system, respectively), rather than splitting off from it. Also considered is a 
related model, for a binary BEC in a single-core trap with the OL, assuming that the two species 
(representing different spin states of the same atom) are coupled by linear interconversion. In that 
case, the symmetry-breaking bifurcations in the AA and RR models switch their character, if the 
inter-species nonlinear interaction becomes stronger than the intra-species nonlinearity. 

PACS numbers: 03.75.Lm, 05.45.Yv, 42.65.Tg, 47.20.Ky 



I. INTRODUCTION AND THE MODEL 



Optical lattices (OLs), i.e., periodic potentials induced by the interference of counterpropagating coher- 
ent laser beams, is a powerful tool that allows one to support various dynamical patterns in Bose-Einstein 
condensates (BECs) I n particular, OLs support gap solitons in BEC with repulsive interaction between 
atoms (the general topic of BEC solitons was reviewed in Ref. Q). Gap solitons in BEC were predicted 
theoretically 0, H[ and then created experimentally in an effectively one-dimensional ( "cigar-shaped" ) trap 
equipped with the OL in the axial direction. In a self-attractive medium, the periodic OL potential makes it 
possible to localize a soliton at a prescribed spot, and supports multi-soliton complexes, as was first demon- 
strated in a model of spatial optical solitons in a planar waveguide with a periodic transverse modulation 
of the refractive index [f| . Similar theoretical results were then reported for solitons in self-attractive BECs 
7|. In particular, in the limit case of a very strong OL, the corresponding Gross-Pitaevskii equation (GPE) 
8 1 reduces to a discrete nonlinear Schrodinger (NLS) equation Q . 

Another topic of great interest to BEC [ltl ] and nonlinear optics (l3j is spontaneous symmetry breaking 
of nonlinear fields trapped in dual-core traps (alias double- well potentials). In particular, the spontaneous 
transition from symmetric solitons in the dual-core trap filled with self-attractive condensate to asymmetric 
solitons is, in the first approximation, tantamount to the formation of asymmetric solitons in dual-core 
nonlinear optical fibers [141 ]. In terms of nonlinear optics, several other dual-core settings have been studied 
theoretically, including ones with quadratic and cubic-quintic [l2j nonlinearities. 

Our objective is to explore symmetric, antisymmetric and asymmetric families of BEC solitons in the 
model of a dual-core trap with cores coupled by linear tunneling and equipped with OLs. We consider 
symmetric systems with attractive or repulsive nonlinearity in each core (to be referred to as AA and RR 
ones, i.e., "attractive-attractive" and "repulsive-repulsive"). In the AA system, the solitons originate in the 
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semi-infinite bandgap of the OL's linear spectrum, while in the RR system there are gap solitons originating 
in finite bandgaps (we will consider the first two gaps, demonstrating that the linear coupling between the 
cores splits them into subgaps). In some cases, soliton families extend into Bloch bands separating the gaps, 
thus becoming embedded solitons [l5j . 

In optics, spontaneous symmetry breaking in two-component gap solitons was studied in models of dual- 
core |16| and triple-core |17| fiber Bragg gratings (solitons in triangular and planar triple-core optical waveg- 
uides without gratings were considered in Refs. [la])- A mixed RA system, with repulsion in one core and 
attraction in the other, will be considered too (the sign of the interaction may be selectively reversed by 
means of the Fcshbach resonance (l9|). 

The basic model that will be dealt with in this work amounts to the following system of linearly coupled 
normalized GPEs for the mean-field BEC wave functions in the two cores, ip(x,t) and <fi(x,t): 

iipt + ipxx + ecos(2x)ip + gi\tp\ 2 ip + kc/> = 

(1) 

i<j>t + 4>xx + ecos(2x + A) 4> + g 2 \4>\ 2 (p + nip = 

where e is the strength of the OL, gx t 2 = ±1 are signs of the nonlincarity k is the linear-coupling coefficient, 
accounted for by the linear tunneling between the cores, and A takes into regard a possible phase shift 
(mismatch) between the OLs in the coupled cores. The above-mentioned symmetry breaking takes place in 
the system with A = 0; a change of the symmetry-breaking bifurcations under the action of the mismatch 
will be considered too. 

In the limit of a very deep OL (e — ► oo), arguments similar to those applied to the single GPE Q suggest 
that the symmetric version of Eqs. (JJ) reduces to a system of two linearly coupled discrete NLS equations. 
In fact, a system of that type was introduced in Ref. [20 ], as a model of a photonic-crystal coupler for optical 
signals. 

The use of Eqs. (JTJ) in the effectively one-dimensional (ID) form implies that both cores are subjected 
to tight transverse confinement, thus allowing one to reduce the underlying three-dimensional GPE to its 
counterparts in ID, as elaborated in several works [21j . Then, the temporal and spatial variables in Eqs. (Q]) 
are related to physical time T and axial coordinate X by 

t = T (n 2 h/md 2 ) , x ee V2irX/d, (2) 

where m is the atomic mass, and d the OL period in physical units. Further, the scaled ID wave functions 
are related to their 3D counterparts, iff and <&, as follows: 



V(X, R,T)\_ to^T I I J 1>(x, t) 1 ( u±m 



$(X,R,T)j-<- y 2\a s \d2\<f>(x,t) j Cxp ( 2h R )' (3) 

where lo±_ and R are the transverse trapping frequency and radial coordinate (around each core), and a s 
the s-wave scattering length of atomic collisions (it is negative if the interaction is attractive). Due to the 
normalizations, the lattice strength is represented by e = Eq/E ycc , where .E rcc = (nh) / (md 2 ) is the lattice 
recoil energy, and E is the depth of the periodic potential in physical units. 

For the experimental realization of the AA and RR systems, respectively, 7 Li [22| and 87 Rb [l[ condensates 
are appropriate. Results for solitons presented below are in the ballpark of n ~ 1 in Eqs. ([1]). With regard to 
Eqs. ((2]), the corresponding tunnel- coupling time, which is t C oupi = 7r/ (2k) in normalized units, translates, for 
d = 1 /jm, into physical coupling time T coup i ~ 10 /is and 100 /is, for lithium and rubidium, respectively. As 
for the number of atoms in the solitons, in the AA system is scales, in the normalized units, between ./V ~ 2 
and N ~ 15, see Fig. (|3|) below. As follows from Eqs. ([2]) and ([3]), this corresponds to the actual number 
of 7 Li atoms per soliton between 10 4 and 10 5 , if experimentally relevant values are assumed, oj± ~ 2ir x 1 
KHz and a s ~ —0.15 nm (in the experiments without the OL, the number of atoms in the solitons was up to 
5, 000 |22|). Similarly, in the RR model the number of 87 Rb atoms in gap solitons is expected to be ~ 1000 



(in the first experiments, it was ~ 250 [l[). 

Equations |T]) may also be interpreted in terms of nonlinear optics, describing spatial distribution of light 
in two parallel planar waveguides, with periodic transverse modulation of the refractive index, i.e., as a 
linear-coupled pair of waveguides considered in Ref. Q. In that case, t is the propagation distance, x is the 
transverse coordinate, while ip and </> are amplitudes of the electromagnetic fields in the planar cores. 
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FIG. 1: Transformation of the linear spectrum with the increase of the lattice strength e. Here and in other 
figures including the spectrum, the shaded areas show Bloch bands, (a) The ordinary (decoupled) Mathicu 
equation, i.e., each equation (j4|) and ([5]) with k = 0; (b) the coupled system, with k = 1 and A = 0. 



Besides the model based on Eqs. ([T]), we will also use a modified one, which additionally includes nonlinear 
coupling between ip and (/>, see Eqs. (|17|) below. It pertains to a single-core trap equipped with the OL, 
that contains a BEC mixture of two different spin states of the same atom [23]. The linear-coupling term 
represents the linear interconversion between the states induced by a spin-flipping electromagnetic wave 
[24l | . While various effects were predicted in the latter setting (25[, the spontaneous symmetry breaking in 
two-component solitons, that we consider in this work, was not studied before. 

The paper is organized as follows. In Section II, we analyze the linear spectrum of the coupled system. 
Methods used for the analysis of soliton solutions are summarized in Section III. Section IV reports results 
for fundamental (even) and twisted (odd) solitons, as well as bound states of solitons, in symmetric systems, 
of both the AA and RR types. Results for asymmetric systems (of the RA type, as well as in systems 
with mismatched OLs) are collected in Section V, including a newly found pseudo-bifurcation, in the system 
with A = 7r. A summary of results obtained in the model which includes nonlinear interaction between the 
coupled components is given in Section VI, and Section VII concludes the paper. 



II. LINEAR SPECTRUM 



Our first objective is to examine how the linear coupling, accounted for by coefficient n in Eqs. (TTJ), alters 
the spectrum of the linear system (first, in the model with aligned lattices, A = 0). We look for solutions to 
the linearized equations as {ip(x, t), (f>(x, <)} = [ct(x) ± (3(x)] e _4Alt , with chemical potential /i. This leads to 
decoupled Mathieu equations (ME) with eigenvalues fj, ± k, 

a" + e cos(2x)a(x) + (fi + n)a(x) = 0, (4) 
13" +ecos{2x)(3(x) + {p- k)(3(x) = 0, (5) 

where the prime stands for d/dx. The ME spectrum gives rise to the well-known bandgap structure. Note 
that, in the model of two linearly coupled Bragg gratings (i.e., a system of four equations for counterpropa- 
gating waves in two cores), the coupling leads to shrinkage or closure of the spectral gap In the present 
case, the effect of the coupling is more complex. Given fj, belongs to a bandgap in the present model if jj, 
falls into one of the gaps in both equations ([?]) and ([5]). Further consideration demonstrates that, depend- 
ing on the values of k and OL strength s, each gap originating from the ME spectrum cither shrinks (or, 
sometimes, completely closes up), similar to the situation in the above-mentioned model of linearly coupled 
Bragg gratings, or splits into pairs of subgaps. An example of the splitting of the first two gaps [Gap {1,2} 
— ► Gaps {(la, lb), (2a, 2b)}] under the action of the linear coupling is displayed in Fig. [TJ 

Unlike Eqs. ([4]) and ([5]), the linearized equations do not decouple with A 7^ in Eqs. ([1]). A typical example 
of the dependence of the spectrum on A is shown in Fig. O With the increase of A, the subgaps (that have 
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FIG. 2: The transformation of the linear spectrum with the increase of mismatch A for e = 8, k = 1. 



appeared as a result of the coupling-induced splitting) shrink. At A close to tt/2, gap 2a disappears, which 
is followed by the disappearance of gap la. At A = w, the entire spectrum qualitatively resembles that of 
the ordinary ME. 



III. ANALYSIS OF SOLITON SOLUTIONS 



A. Symmetric, antisymmetric, and asymmetric solitons 



Stationary solutions to Eqs. |T|) are looked as {ip, 4>} = {u(x),v(x)} e ltit , with real functions u and v 
obeying the equations 



fxu + u" + e cos(2x)u + giu 3 + nv = 0, 
[iv + v" + e cos(2a; + A)v + g2V 3 + ku = 



3^„._n ( 6 ) 



First, we explore the symmetries of the system in the AA (gi = g2 = +1) and RR (gi = 92 = — 1) cases 
with aligned OLs, A = 0. To this end, we designate uq(x; fi) the solution of the single-component GPE, i.e. 
either of two equations (J6j) with k = A = 0, corresponding to chemical potential /1. Then, Eqs. ([6]) give rise 
to symmetric and antisymmetric solitons, 

u(x; jj) — v{x\fi) = Uq(x; (J, + k), (7) 



u(x;n) = -v(x;fj,)=uo(x;fJ,-K). (8) 

If the linear coupling splits the bandgap into two subgaps as outlined above, the shifts of the chemical 
potential as per Eqs. (J7]) and l[5]). fi — > fi ± k, displace the symmetric and antisymmetric soliton to the 
lower and higher subgap, respectively. In particular, in the coupled AA system, symmetric solitons, which 
are generated by their counterparts belonging to the semi-infinite gap in the single-component GPE with 
attraction according to Eq. l[7|). remain in the semi-infinite gap, while the corresponding anti-symmetric 
solitons may either stay in the semi-infinite gap, or move to subgap la (in the case shown in Fig. [1]), or 
even end up inside the Bloch band separating the semi-infinite gap and subgap la. In the latter case, the 
antisymmetric soliton becomes embedded (into the continuous spectrum), as defined in Refs. [l5j . 

The central issue to be considered below is a possibility of spontaneous symmetry breaking, i.e. finding a 
bifurcation point beyond which asymmetric solitons emerge. To quantify the symmetry breaking, we define 
norms in the two cores, and the asymmetry ratio, O, as 

N u , v = f~{u*(x)y(x)}dx, 9 5 ^g (9) 

[Eqs. {1} conserve only the total norm, N = N u + N v ]. The bifurcation was found from a numerical solution 
of Eqs. ([6]). The results are reported in the next section for the models of the AA, RR, and RA types. 
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B. Soliton stability 



To tackle the stability problem, perturbed solutions are taken as 

^(ar.i)} = {[u(x ]f i)+^(x)e tnt ] ,[v(x;ijl) + 60^'] } e^\ 

(M)} = {[fi*(x;/i) + 6(x)c^] ,[V(x;v) K"*. 



(10) 



where /x), //)} e~ lflt , with real functions u(x; /i) are £(x;/x), are stationary solutions of Eqs. 

Perturbations ^ 3 of fields ip,4>, and their complex conjugates, £2,4, are formally treated as independent 
functions. Substituting expressions (|10p in Eqs. (|TJ) and linearizing, we arrive at an eigenvalue problem, 



£(u(x;/x),/x,#i,0) kct 3 

k<t 3 C(v(x',fi),fi,g2,A) 





"62" 


= 


62" 




.6* 


&4. 



where 62 = (6,6) ,64 = (6,6) 
£(/(x),/x,<7,A) = 



the linear-stability operator for the single-component GPE is 
d 2 



dx 2 



e cos(2x + A) + 2gf 2 {x) + /x 



0-3 + (x)tT2, 



(11) 



(12) 



and (73 and 02 are the Pauli matrices. In the single-component GPE, families of fundamental solitons have 
been shown to be stable [2(|, i.e., all eigenvalues of the corresponding operator C (uq{x; /x), fx, g, 0) are real. 

Equation (JTTJ) can be simplified for symmetric and anti-symmetric solutions. Substituting, respectively, 
expressions (J7J and ((SJ) (with A = and 51 = 92 = ff), we arrive at a reduced eigenvalue problem for the 
symmetric soliton, 



C (u {x; fj, + K),fj, + K,g,0)^ + = 776-, 
C(uo(x;(i + k),(i- K,g,Q)£- = 

(6r = 62 i 64)1 an< l its counterpart for the anti-symmetric one, 

C(u (x;fj,- K),(j, + K,g,Q)£ + = ?^+, 
xC(m (x;/x- k),/x- k,5,0)6- = »7C— 



(13) 
(14) 



(15) 
(16) 



Since Eqs. (|T5)| and (|T5j) are tantamount to the linear-stability problem in the single-component GPE, the 
stability of the soliton in the latter equation is a necessary condition for the full stability of symmetric and 
antisymmetric solitons in the coupled system. To identify the full stability conditions for the symmetric and 
antisymmetric solitons, we solved the additional eigenvalue problems, i.e., respectively, Eqs. (|14p and (|15[) . 
The stability of asymmetric solitons was inferred from a numerical solution of the full eigenvalue problem, 
based on Eq. (fTTj) . The solution is stable if all respective eigenvalues 77 are real. 



IV. RESULTS (SYMMETRIC SYSTEMS) 
A. The attraction-attraction model 

Families of numerically found stationary soliton solutions of Eqs. © for the AA system (g\ = #2 = +1) 
with zero mismatch (A = 0) are displayed in Fig. [3] for k = 1, in the cases of weak (e = 1) and strong 
(e = 8) OLs. We observe that a branch of asymmetric solitons bifurcates from the symmetric one, while 
the antisymmetric solutions do not give rise to any bifurcation. We also notice that the antisymmetric 
branch extends through the Bloch band separating the semi-infinite gap and the first finite bandgap, la; as 
said above, the solitons are embedded ones inside the band. Panels (b) and (d) in the figure clearly show 
that the symmetry-breaking bifurcation is supercritical (the bifurcating branch immediately goes forward 
in — /x, without turning backward, the latter being a characteristic feature of the subcritical bifurcation - 
in particular, in the system of linearly coupled NLS equations of the AA type without the lattice potential 

14 1). The present situation is, generally, similar to that in the system of linearly coupled Bragg gratings 

1 61 ] (where the bifurcation is supercritical too). 
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FIG. 3: (Color online) Soliton families in the attraction-attraction model with n = 1. In this figure and 
below, the results are displayed for the dual-core model with aligned lattices (A = 0), unless a nonzero 
value of A is indicated. Presented are cases of relatively weak, e — 1 (a.b), and strong, e — 8 (c,d), lattices. 
Here and below, stable and unstable branches are shown by solid and dashed lines, with labels "s", "an" 
and "a" referring to symmetric, antisymmetric, and asymmetric solutions, respectively. (a,c): The total 
norm of the soliton versus its chemical potential, fj,. (b,d): The asymmetry ratio [defined as per Eq. ©] 
versus /i, for families of asymmetric solitons; the vertical line labeled ^bif marks the bifurcation point. Note 
that the branch of anti-symmetric solutions is completely unstable in the weak lattice (a,b), and completely 

stable in the strong one (b,d). 



The change of the bifurcation with variation of coupling constant n is illustrated by Fig. 2] In particular, 
the critical value of the coupling constant, Kbifj at which the bifurcation occurs is shown in panel (b) as 
a function of the soliton's total norm [note that in Fig. GJa), which pertains to n = 1, the bifurcation in 
the system with e = 1 occurs at N f=a 4, in compliance with Fig. Bib)]. Naturally, «bif grows with N, as 
the bifurcation is a result of the competition between the nonlinear and linear properties of the system, 
accounted for by N and k, respectively. On the other hand, the dependence of the bifurcation point on 
the lattice's depth e is quite weak, the stronger lattice making the region of the existence of asymmetric 
solitons somewhat larger (which can be easily understood too, as the lattice tends to pin and thus stabilize 
any localized pattern). 

The above results bear some implications for stability of the solutions, as, in models of the NLS type with 
self- focusing, a necessary stability condition is given by the Vakhitov-Kolokolov (VK) criterion [28|. In the 
present notation, it is dN/d^i < 0. As seen in Fig. [3J all solution branches satisfy this condition (however, 
as shown below, not all of them are stable, as the VK criterion is not sufficient for the stability). 
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FIG. 4: (Color online) Properties of asymmetric solitons in the attraction-attraction system, (a) The 
asymmetry ratio [defined as per Eq. ([9])] as a function of coupling coefficient k at e — 5 and different values 
of the soliton's norm, (b) The value of n at the bifurcation point versus the soliton's norm at different 

values of the lattice strength, e. 
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FIG. 5: (Color online) The largest instability growth rates for symmetric and anti-symmetric solitons in 
the attraction-attraction system, with e = 1, k = 1, are shown by the continuous and dashed curves, 

respectively. 



To investigate the stability of the symmetric and antisymmetric soliton families presented in Fig. [3] in 
an accurate form, we solved the reduced eigenvalue problem, based on Eq. (Ti"4|) or (p~5|) . respectively, to 
identify the eigenvalue with the largest imaginary part, Im{ry}, i.e., the fastest growing unstable mode. 
Figure [5] shows the result for a relatively weak OL, with e = 1. We conclude that the symmetric soliton is 
stable up to the bifurcation point, where the branch of the asymmetric solitons emerges. Beyond this point, 
the symmetric solitons are unstable. On the other hand, the branch of the anti-symmetric solitons, which 



undergoes no bifurcation in the AA system, is unstable in weak OLs, see Figs. 3(a) and [5l but becomes 



stable in a stronger lattice, see Fig. 3(c) In fact, strong OLs give rise to bistability: the antisymmetric 
soliton coexists, as a stable solution, with either symmetric or asymmetric one. Finally, numerical solution 
of the full eigenvalue problem, based on Eq. (jllj) . for the asymmetric solitons demonstrates that they are 
stable whenever they exist. 

The predictions of the linear stability analysis have been confirmed by direct numerical simulations of the 
underlying equations ([T]). In particular, Fig. [B] shows that unstable symmetric solitons evolve into stable 
asymmetric ones. It has also been checked that antisymmetric solitons in sufficiently strong OLs are stable 



FIG. 6: (Color online) Perturbation-induced transformation of an unstable symmetric soliton into a stable 
asymmetric one in the attraction- attraction model with a relatively weak lattice (e = 1, k = 1). The initial 
norm is N = 5. 6. Here and in Fig. 1141 below, the two components of the solution are juxtaposed in the same 

panel. 
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FIG. 7: (Color online) The attraction-attraction model: critical values of the coupling coefficient, K cr , 
versus the lattice strength e, for a fixed norm, N = 10. The upper curve (squares) is, in fact, the 
bifurcation line, i.e., the stability border of the symmetric solitons (which are stable above the border), and 

existence border of stable asymmetric solitons (they exist beneath the border). The lower curve 
(rhombuses) is the stability border of antisymmetric solitons, which are stable beneath it, thus making the 

lower area a bistability region. 



indeed, while, in weak lattices, they are completely destroyed by growing perturbations. 

Results of the analysis are summarized in Fig. [3 in the form of a diagram showing the stability regions 
of the symmetric, antisymmetric, and asymmetric solitons in the (e, k) plane. To this end, we define, in 
addition to the above-mentioned bifurcation value, Kbifi which plays the role of the instability border of 
symmetric solitons (they are unstable, for given e and N, at k < Kbif, and stable at k > Kbif) and existence 
border of stable asymmetric solitons (they exist at k < Kbit), a critical value which separates stable and 
unstable antisymmetric solitons: they are stable at k < K cr , and unstable at k > n cr . 

The stability analysis demonstrates that the branch of the antisymmetric solitons changes its character 
from unstable to stable as a whole, with the increase of e at fixed k [see Figs. [3ja) and (c)], or decrease of 
K at fixed e. In other words, a situation has not been found when a part of the branch of antisymmetric 
solitons would be stable, while its other part is unstable. Note that the segment of the branch extending 
through the Bloch band, where the antisymmetric solitons are embedded ones, is stable too. 
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FIG. 9: (Color online) Continuation of Fig. [5] to larger values of the soliton's norm. 



B. The repulsion-repulsion system 

In the single-component GPE with repulsion, solitons may only be found in finite bandgaps, therefore 
they are called gap solitons [H, I n the coupled system with repulsion in both cores, symmetric and 
antisymmetric states may be composed of single-component gap solitons as per Eqs. ([7]) and (jSJ). Families 
of the solitons in the RR system, generated from the gap solitons belonging to the first and second finite 
bandgaps in the single-component GPE, are shown in Figs. [5] and [51 along with families of asymmetric 
solitons bifurcating from them. In particular, the symmetric branch in subgaps 2a and 2b in Fig. [9] (which 
extends across the Bloch band separating the two subgaps), and antisymmetric branch in subgap 2b (it 
extends into the adjacent Bloch band) may be regarded as continuations of the, respectively, symmetric 
branch in subgap lb, and antisymmetric one in subgap 2a, which are displayed in Fig. [5J On the other 
hand, the origin of the stable branch of asymmetric solitons in subgap 2b in Fig. [8] is not clear, as numerical 
problems impede to continue it in the direction of decreasing fi and N. 

On the contrary to the situation in the AA model [cf. Figs. [3] (a) and (c)], in the present case the 
symmetric branch does not undergo any bifurcation, while asymmetric solitons bifurcate supcrcritically 
from the antisymmetric branch. Dependences of asymmetry ratio O on k, and of the bifurcation value, Kbif , 
oniVs N u + N v for asymmetric solitons in the RR system are quite similar to those in its AA counterpart, 
which were shown above in Fig. 2] (therefore, the dependences for the RR system are not displayed here). 
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FIG. 10: (Color online) An asymmetric soliton in the repulsive-repulsive system, with e = 8, K = 1, 

N u = 19, N v = 0.9, n = 5.39. 




FIG. 11: (Color online) The instability growth rate for antisymmetric solitons in the first two (sub)gaps in 

the repulsion- repulsion model with k — 1 and e = 8. 



A typical shape of the asymmetric soliton is presented in Fig. [TO) It is worthy to note that both symmetric 
and antisymmetric soliton branches, but not asymmetric ones, may become embedded, crossing Bloch bands 
which separate the subgaps. 

Similar to what was reported above for the AA model, the reduced eigenvalue problem based on Eqs. (|14[) 
and ()15j) was solved numerically to determine the stability of symmetric and antisymmetric gap solutions 
in the RR system (the VK criterion is irrelevant for self-defocusing models). The analysis demonstrates 
that the families of symmetric solitons presented in Figs. [5] and [5]) are completely stable if the OL is strong 
enough, but symmetric solitons may be unstable in a weak lattice (for instance, at e — 1). This property 
is reminiscent of the AA system, where the antisymmetric solitons are unstable at e = 1, and stable at 
large e, see Figs. [31 E] and [7] Further, antisymmetric solitons are stable before the (anti)symmetry-breaking 
bifurcation, i.e., at fx < Hha, and, quite naturally, unstable at /i > /ibif (even if the antisymmetric branch 
does not seem to be continuously connected to the bifurcation point, see Fig. [9]). In the latter case, the 
instability growth rate for the antisymmetric solutions is shown in Fig. 1111 

Asymmetric solitons bifurcating from the antisymmetric ones are stable whenever they exist, which was 
checked by solving the full eigenvalue problem for them, based on Eqs. (jlip . As well as in the case of the 
AA system, these results imply bistability, provided that the OL is strong enough. Indeed, the symmetric 
solitons are then stable, and, along with them, either antisymmetric solitons or asymmetric ones (below or 
above the bifurcation, respectively) are stable too. 

The predicted stability and instability of symmetric, antisymmetric, and asymmetric solitons has been 



11 




|X x 

(a) (b) 

FIG. 12: (Color online) (a) Families of twisted (odd) solitons in the attraction-attraction system, with 
e = 8 and k = 1. Symmetric twisted solitons are stable up to the bifurcation point, beyond which they are 
unstable, while the emerging twisted asymmetric solitons are stable. Antisymmetric twisted solitons are 
stable everywhere (actually, because the lattice is strong enough), cf. Fig. [3Jc). (b) An example of an 
asymmetric twisted soliton with N u = 5.3, N v = 0.7, fi ~ —0.66. In this figure and in some figures below, 
shaded vertical stripes depict the lattice potential which supports the solitons. 



verified by direct simulations too. In particular, the antisymmetric solitons, if unstable, tend to rearrange 
themselves into their stable asymmetric counterparts. 



C. Twisted solitons and bound states 



Twisted (alias odd) solitons, which feature two out-of-phase amplitude peaks in one period of the under- 
lying lattice, were found and shown to be stable in the single-component GPE with attractive nonlinearity 



26 1 (unlike formally similar subfundamental gap solitons in the model with self-repulsion, which are unstable 

We have studied spontaneous symmetry breaking of twisted solitons in the AA system. Families of twisted 
solitons in this system are shown in Fig. I12f a). The behavior is generally similar to that of the fundamental 
solitons in the same (AA) model: a stable asymmetric branch bifurcates from the symmetric one, leaving it 
unstable; however, it is worthy to note that, in the case shown in Fig. [T27 a). the bifurcation occurs in the 
finite (sub)gap, lb, while the similar bifurcation of fundamental solitons was observed, in Figs. [3]Ja) and 
(c), in the semi-infinite gap. The entire antisymmetric branch of the twisted solitons is stable, provided that 
the lattice is strong enough, or unstable otherwise. A typical example of the asymmetric twisted soliton is 
shown in Fig. [H^b). 

Bound states of 7r-out-of-phase fundamental solitons, which also feature the odd parity, but include peaks 
separated by an empty lattice site, have been found too, in the AA and RR systems alike (in the single- 
component GPE with the repulsive nonlinearity, the presence of an empty site between the peaks is a 
necessary condition for the existence of a stable bound state of two out-of-phase fundamental gap solitons 
{27ll2"9l|). Families of the bound states of this type are shown in Fig. [T3l being quite similar to fundamental- 
soliton families, cf. Figs. [3]and[5J A difference is observed in the evolution of symmetric and anti-symmetric 
solitons destabilized by the bifurcation: rather than rearranging into asymmetric solitons, they give rise to 
persistent breathers, as shown in Fig. 1141 



V. ASYMMETRIC SYSTEMS 



As said in Introduction, the two-core system can be made asymmetric by either assuming that the signs 
of the scattering lengths arc opposite in the parallel-coupled traps, or by admitting a mismatch between the 
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FIG. 13: (Color online) Families of twisted bound states, for e = 8, k = 1. Insets show typical 
bound-soliton profiles. Panels (a) and (b) pertain to the attraction-attraction and repulsion-repulsion 

system, respectively. 





(a) 



(b) 



FIG. 14: (Color online) Evolution of unstable twisted (odd) bound states of fundamental solitons in the 
system with e = 8, k = 1. (a) A bound state of symmetric solitons in the attraction-attraction system; (b) 
a bound state of antisymmetric solitons in the repulsion-repulsion system. In either case, formation of a 

long-lived breather is observed. 



two OLs. Both cases are considered in this section. 



A. The repulsion-attraction system with aligned lattices 

In the RA system (we set g\ = —1, cj2 = 1, assuming opposite signs but equal absolute values of the 
nonlinear coefficients in the components), numerical solution of Eqs. ^ reveals two soliton families, viz., 
ones with a dominant repulsive component (N u > N v ) residing in subgap lb, or with a dominant attractive 
component (N u < N v ) in the semi-infinite bandgap (solitons in two-core RA systems, but without the 
lattice, were considered in Refs. [I(|). Generic examples of such families are presented in Fig. [1^1 and 
typical examples of the respective soliton shapes are displayed in Fig. [16] 

The computation of the eigenvalues for perturbation modes demonstrates that both families are completely 
stable. Note that the stability of the attraction-dominated solutions, with N u < N v , is also suggested by 
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FIG. 15: (Color online) Solitons families in the repulsion-attraction model (—gi = .92 = 1) with e = 8 and 

K= 1. 
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FIG. 16: (Color online) Examples of solitons in the repulsion-attraction model, with —g\ = g% = 1, and 
e = 8, k = 2. (a) iV u = 5.1, i\T„ = 0.9, fj, w -0.105; (b) A u = 0.6, iV„ = 5.4, /z w -9.58. In either case, the 

total norm is TV = 6. 



the VK criterion, if applied to dependence N(/j,) in Fig. I15f a) (for the repulsion-dominated solutions, with 
N u > N v , the VK criterion is irrelevant, as mentioned above). 



B. Bifurcations and pseudo-bifurcations in systems with mismatched lattices 

To consider effects of the mismatch between the OLs in the two cores, we will here focus on Eqs. (JTJ) with 
A = 7r/2 and A = 7r (the latter value corresponding to the largest mismatch, similar to a system of parallel 
Bragg gratings with the maximum mismatch, that was recently considered in Ref. [3l[). For all types of 
the nonlinear interactions in the cores, at k = the solitons in the decoupled shifted lattices were taken 
as u(x) = Uo(x), v(x) = 0, where uq(x) is the corresponding soliton in the single-component GPE. At this 
point, the asymmetry ratio defined in Eq. ([9]) is at its maximum, = 1. As the linear-coupling coefficient, 
k, increases, the solution becomes less asymmetric. 

Figure [T7] displays soliton families found in the AA system. As in Ref. [12], which treated a mismatched 
system of parallel-coupled Bragg gratings, we distinguish between asymmetric and quasi-symmetric (QS) 
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solitons [a branch of quasi-antisymmetric (QAS) solutions was found too but is not shown, as it is completely 
unstable, cf. Fig. [3J. QS solitons have equal norms in both components, i.e., O = 0, and feature similar, 
although not identical, profiles of the components. The bifurcation chart for the AA system with A = 7r/2 
is very similar to the one for the same system with aligned OLs (A = 0), the symmetric solutions being 
replaced by their QS counterparts. In particular, the asymmetric branch bifurcates from the QS one. 

At first glance, it may seem that, in the AA system with A = n, the branch of asymmetric solitons 
also bifurcates from the QS one. However, this is not the case. A blow-up (inset in Fig. [T7b) shows a 
pseudo-bifurcation, which means that asymmetric and QS branches gradually approach each other to a point 
where they seem indistinguishable, but they never merge, therefore the QS solutions always remain (strictly 
speaking) unstable, while the asymmetric ones are always stable (the stability is discussed in more detail 
below). In Fig |18i we display an example of comparison between shapes of asymmetric and QS solitons with 
equal norms, N = 3.5, when they are very close to each other, the asymmetric one having m 10 -4 (its QS 
counterpart has 8 exactly equal to zero). Note that the centers of the two components in the QS soliton are 
slightly separated, while in the asymmetric soliton they exactly coincide. 

Figure [TO] additionally illustrates the difference between the true quasi-symmetry-breaking bifurcation in 
the AA system with A = tt/2, and the pseudo-bifurcation in the system with A = ir. In the latter case, the 
branches closely approach but never merge, cf. Fig. [T7] 

Typical shapes of strongly asymmetric solitons in the mismatched AA model are displayed in Fig. 1201 
Naturally, the solitons are more asymmetric at A = it . 

Families of solitons found in the mismatched RR system are presented in Fig. [2U The QS branch is 
called this way only in the sense that peak amplitudes of the two components have the same sign. In fact, 
profiles of the components in the QS soliton may be very different, as insets show in Fig. 1211 The quasi- 
antisymmetry-breaking bifurcation in the RR system and respective pseudo-bifurcation are illustrated by 
Fig. [221 and a set of typical profiles of the solitons is displayed in Fig. [221 

Stability of the solitons in the mismatched AA and RR systems was studied in direct simulations of Eqs. 
([T]) over a wide range of values of N and k. In either case (as in the aligned system), the asymmetric solitons 
were found to be stable whenever they exist, and the branch which gives rise to the bifurcation is stable 
before the bifurcation and unstable afterwards. Further, QAS and QS solitons in the AA and RR system, 
respectively, are unstable at small e and stable at larger e. When QS solitons are unstable in the AA system, 
they evolve into their stable asymmetric counterparts. Similarly, in the misaligned RR system, QAS solitons 
which were destabilized by the bifurcation (or which are unstable due to the pseudo-bifurcation) transform 
themselves into stable asymmetric solutions. 

Soliton families in the mismatched system of the RA type were investigated too. It was found that the 
corresponding picture is quite similar to that presented for the RA system without mismatch in Fig. 1151 
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FIG. 18: (Color online) Comparison of soliton profiles in the attraction-attraction model with A = ir and 
e = 1, k = 1. (a) A quasi-symmetric soliton, with = 0; (b) a slightly asymmetric soliton, with w 10 -4 . 
Both solitons have equal norms, TV = 3.5. Here and below, the set of vertical shaded stripes represents the 

periodic lattice potential in the first core. 
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FIG. 19: (Color online) The symmetry ratio for solitons in the mismatched attraction-attraction model 
with e = 1 versus the coupling coefficient, k: (a) A = n/2; (b) A = n (similar dependences for the aligned 

system, with A = 0, are displayed in Fig. |4j. 



VI. BINARY MIXTURES WITH LINEAR COUPLING 



In this section we consider the model of a binary BEC in the single-core trap equipped with an OL. As 
explained in Introduction, the two components of the mixture (ip and </>) represent two different spin states 
of the same atom, and the linear coupling between them is induced by a spin-flipping electromagnetic wave. 
The corresponding system of the normalized GPEs for wave functions ip and <f> takes the form 

iik + i>xx + £ cos(2a:)V> + (gi\ip\ 2 + 512H 2 ) tp + K(j) = 0, , . 

i<pt + <Pxx+£cos{2x)(j)+{g2\<iy\ 2 +gi2\i>\ 2 )<P + ^ = U K ' 

(in the single trap, there may be no mismatch between the lattice-potential terms in the two equations). As 
the interaction between atoms may be both attractive and repulsive, in this section we consider two basic 
cases: when all the three nonlinear terms are attractive, i.e., 171 = g% = 1, .912 > (AAA), and when they all 
are repulsive, with g\ = 92 = — 1, .912 < (RRR). 
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FIG. 20: (Color online) Typical asymmetric solitons in the mismatched attraction-attraction system with 
e = 1: (a) A = tt/2, N = 2.5, 6 = 0.51, fj, = -0.83, K = 0.4; (b) A = tt, N = 2.5, 9 = 0.63, fx = -0.85, 

K = 0.4. 
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FIG. 21: (Color online) Soliton families in the misaligned repulsion-repulsion system, with e — 8. Labels 
"a" , "qs" , and "qas" refer to asymmetric, quasi-symmetric and quasi-antisymmetric branches, respectively. 
Insets show soliton profiles at points marked by diamonds, (a) A = 7r/2; (b) A = 7r. 

Bifurcation diagrams for the AAA and RRR systems are shown in Fig. [2H In the former case (AAA), 
stable asymmetric solitons bifurcate, as before, from symmetric ones if < g\2 < 1. However, if the 
inter-species nonlinear interactions dominate, i.e., for gi2 > 1, the character of the bifurcation changes, 
and asymmetric solitons bifurcate from antisymmetric ones. In the RRR system, a similar transition is 
observed: the asymmetric solitons are generated, as above, by the bifurcation from antisymmetric solitons if 
< —.9i2 < 1, but, for —,gi2 > 1, antisymmetric solitons bifurcate from symmetric ones. In the Manakov's 
case, | <7i 2 1 = 1 [H|, no bifurcation was found. Note that linearly coupled equations (jTTJ) without the lattice 
potential (e = 0), but with the linear coupling present, k ^ 0, are equivalent to the Manakov's integrable 
system with n = [34| , therefore the absence of the bifurcation in this case is not surprising. 

Direct simulations of the AAA system demonstrate that the asymmetric solitons, whenever they exist, are 
stable, and, as usual, the bifurcation destabilizes the solitons from which the asymmetric ones emerge. In 
the RRR system, the stability situation is more complex. While asymmetric solitons bifurcating from the 
symmetric ones in the RRR system, i.e., in the case of | <7i2 1 > 1, are always stable, the asymmetric solitons 
bifurcating, in the same system, from the antisymmetric branch (at < | <7i2 1 < 1) may be both stable and 
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FIG. 23: (Color online) Profiles of asymmetric solitons in the mismatched repulsion-repulsion model with 
e = 8. (a) A = tt/2, N = 2, 9 = 0, fx = -2.81, k = 1.5; (b) A = tt/2, N = 10, 9 = 0.31, /i = 0.23, k = 2.5; 
(c) A = 7T, = 2, 9 = 0.92, /j, = -2.96, k = 1; (d) A = tt, JV = 10, 6 = 0.13, /i = -5.57, re = 8. 
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FIG. 24: The bifurcation value of the coupling coefficient, Kbif > versus the relative strength of the 
interspecies nonlinear interaction, g±2, in the model of the binary mixture loaded into the optical lattice 
with strength e = 8 [Eqs. (|17|) ]. The vertical dashed line divides regions in which asymmetric solitons 
bifurcate from symmetric and antisymmetric ones, (a) The model with attractive nonlincarities (for total 
norm TV = 10). In this system, the asymmetric solitons are always stable, (b) The model with repulsive 
nonlincarities (for TV = 4). In this case, all asymmetric solitons bifurcating from the symmetric ones are 
stable, while the solitons bifurcating from the antisymmetric branch are unstable if \g\2\ is close to 1 (but 

they are stable if \g^\ is small). 



unstable. We did not try to identify a border between the stable and unstable asymmetric solitons in the 
latter case in an accurate form; however, the asymmetric solitons in the RRR system are definitely stable 
for |c/i2| — > (which complies with the results reported above for the RR system with g\2 = 0), and they are 
unstable for \gn\ close to (but smaller than) 1. In that case, direct simulations demonstrate that unstable 
asymmetric solitons transform into persistent breathers (not shown here). We did not systematically study 
the stability of the nonbifurcating families, i.e., antisymmetric and symmetric ones in the AAA and RRR 
systems, respectively. 



VII. CONCLUSION 



In this work, we have studied families of soliton states in a symmetric set of two effectively onc-dimcnsional 
traps (cores) equipped with OLs (optical lattices), filled with self-attractive or self-repulsive BEC, and 
coupled by linear tunneling (the same model may be also realized in terms of spatial solitons in two parallel 
planar optical waveguides, which carry lattices in the form of a transverse modulation of the refractive 
index). Asymmetric systems, with a phase shift A between the two OLs, as well as with opposite signs 
of the nonlincarity in the BEC confined in the different traps, were considered too. In the spectrum of 
the dual-core system, the linear coupling between the cores splits finite bandgaps induced by the single OL 
into subgaps, or partly closes them. Influence of mismatch A on the spectrum was considered too, with a 
conclusion that large A tends to make the spectrum qualitatively similar to that in the single-core model. 
In the full nonlinear model, families of symmetric and antisymmetric solitons, as well as asymmetric solitons 
generated by symmetry-breaking bifurcations, have been constructed. This was done in the case of attraction 
or repulsion in both cores (AA and RR systems), as well as in the mixed (RA) system. Deformation of the 
soliton families in the AA and RR system under the action of mismatch A between the two OLs was 
investigated too. 

In the AA and RR systems, asymmetric solitons bifurcate supercritically from symmetric and antisym- 
metric ones, respectively. In the former case, symmetric and asymmetric fundamental solitons were found 
only in the semi-infinite gap, while antisymmetric solitons were also discovered in the first finite subgap, 
as well in the Bloch band separating it from the semi-infinite band (in the latter case, the antisymmetric 
solitons have the character of embedded ones). In the RR system, solitons are located in finite bandgaps, 
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and, additionally, antisymmetric and symmetric solitons were found as embedded ones inside Bloch bands. 
In both systems (AA and RR), the asymmetric solitons are always stable, while the branch which gives 
rise to them is stable before the bifurcation, and unstable afterwards. The nonbifurcating branches, i.e., 
antisymmetric and symmetric solutions in the AA and RR systems, respectively, are unstable or stable as a 
whole in weak and strong OL, respectively. The corresponding stability region was identified in a complete 
form for the AA system, in the parameter plane of coupling coefficient k and OL strength e. In the systems 
of both AA and RR types, bistability regions were found, in which the stable asymmetric solitons coexist 
with cither antisymmetric or symmetric ones, respectively. If the symmetric soliton in the AA system, or the 
antisymmetric soliton in the RR system was destabilized by the bifurcation, its evolution tends to transform 
it into a stable asymmetric soliton, in either system. In the RA system, two soliton families were found, 
dominated by either the self-attractive or repulsive component, both being entirely stable. 

In the AA system, we have also studied families of twisted (odd) solitons, featuring two out-of-phase peaks 
in a cell of the OL. The behavior of these families is, generally, similar to that of the fundamental solitons, 
including the bifurcation of the symmetric branch into an asymmetric one, and bistability; however, unlike 
the situation with the fundamental solitons, the symmetry-breaking bifurcation of the twisted solitons occurs 
in a finite bandgap. In addition, 7r-out-of-phasc bound states of fundamental solitons were constructed in 
both the AA and RR systems. Their behavior too is quite similar to that of the respective fundamental 
solitons, with a difference that, past the bifurcation, an unstable bound state tends to arrange itself into a 
persistent breather (rather than into a stationary bound state of the asymmetric type). 

The study of the systems with the mismatch, A, between the OLs in the two cores demonstrates that former 
symmetric and antisymmetric solitons turn into QS (quasi-symmetric) and QAS (quasi-antisymmetric) ones. 
For A = 7r/2, the deformation is mild, and the behavior of various soliton families remains qualitatively the 
same as in the system with A = 0; in particular, well-defined bifurcations which generate asymmetric solitons 
from QS or QAS ones, in the mismatched AA and RR system, respectively, occur despite the presence of 
the mismatch (in a system of mismatched parallel-coupled Bragg gratings, similar observations were recently 
made in Ref. [13|). However, when the mismatch attains its maximum, A = n, the situation becomes 
different: instead of the bifurcations, both the AA and RR systems demonstrate pseudo-bifurcations, when 
the branch of asymmetric solutions does not really bifurcate from that of QS or QAS solutions, but rather 
goes very close and only asymptotically merges into it. Thus, the asymmetric branch always exists and 
is stable in the latter case, while its QS or QAS counterpart is always unstable (although the instability 
is extremely weak when it approaches the stable asymmetric branch). To the best of our knowledge, the 
present work reports the first example of the pseudo-bifurcation in a mismatched binary system (in a system 
of mismatched Bragg gratings, no such effect was observed (32|V 

Finally, we have considered a model of the single-core trap, equipped with the OL and filled with a binary 
mixture, assuming both the nonlinear interaction between the two species, and linear coupling between 
them, induced by a resonant spin-flipping field, if the species correspond to different spin states of the same 
atom. In this case, we have found that the bifurcations generating asymmetric solitons from symmetric 
and antisymmetric ones are qualitatively similar to their counterparts in the model without the nonlinear 
interaction between the species, if this interaction is weaker than the intra-species nonlinearity; otherwise, 
the bifurcations switch their character, so that the asymmetric solitons emerge from antisymmetric and 
symmetric ones in the AA and RR systems, respectively. 

A very natural extension would be to consider a two-dimensional version of the models introduced in this 
work, which may also directly apply to BEC in a dual-core pancake-shaped trap, or to a binary mixture 
in a single-core nearly flat trap. The work in this direction is currently in progress, and results will be 
reported elsewhere. It may also be possible to consider the spontaneous symmetry breaking in a degenerate 
fcrmion gas filling a dual-core trap, using a relatively simple description of the fermion gas based on the 
Thomas-Fermi (alias mean-field-hydrodynamic) approximation, which was developed, for various settings, 
inRefs. d|. 
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